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1. INTRODUCTION 


1.1 THE FINITE ELEMENT METHOD 

The finite element method, in general, is an approximate method to 
solve differential equations. Using variational calculus the differential 
equation under consideration is posed as a functional. The resulting 
functional depends upon the unknowns and their derivatives with respect to 
the spatial coordinates x,y and z and possibly the time, t. In structural 
problems the functional represents a meaningful quantity, namely, the 
potential energy. However, in general, the functional may not have any 
physical interpretation. Minimizing the functional with respect to the 
unknowns is equivalent to solving the differential equation. The functional 
is minimized by setting its first variation to zero. In structural problems this 
corresponds to the well known concept of minimization of the potential 
energy. The result of the minimization is a set of algebraic equations 

(1.1) [K]{u} = {f> 

where [ K ] is the matrix of coefficients of the unknowns and is known as the 
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"stiffness matrix", 

{ u } is an array of unknowns and 
{ f } is an array of forcing functions. 

The equations are then solved for u. In a typical application the 
domain under consideration is modeled by dividing it into elements. An 
interpolation function, or shape function, is set for the elements to interpolate 
values of the unknown at any point inside the element in terms of its values 
at the nodes. This interpolation function is used in the functional which when 
minimized as described above, yields the stiffness matrix. The reader is 
directed to several excellent books on finite element method by Zienkiewics 
[1,2]*, Segerlind [3], Reddy [4], and Huston [5]. The point to note for this 
report is the important role of the minimization process involved in the finite 
element methods. 


12 MESH GENERATION 

Each element in the finite element model is addressed by its 


•Numbers in brackets indicate references listed at the end. 
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number. Also each node is addressed by its number. The inter-connectivity 
of the elements is determined by the common nodes shared by the 
elements. In a model with few elements and nodes, the user can manually 
divide the domain, number each element and node, and keep track of the 
element connectivity. However, in models with many nodes and elements, 
the effort required to divide the domain into elements and attend to 
connectivity is great. It then becomes difficult to accomplish this task 
without committing errors. However, there are several finite element pre- 
processors which do this job automatically once the geometry is defined. 
Users can then devote more time to interpreting results. Shephard [6] 
has reviewed the current trends in mesh generation. Although there are 
several ways to generate meshes, these methods fall into two broad 
categories: 


1.2.1 MAPPING TECHNIQUES 

This type of mesh generation is best suited when the geometry is 
simple - as in the case of a rectangle or a cuboid. Typically the user 
needs to choose the number of elements on each of the edges that 
defines the geometry and the element concentration along the edges. The 
software then generates the mesh simply by joining nodes on the opposite 
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edges. NASTRAN MSGMESHf , GIFTSd , SUPERSAP* and 
SUPARTABt (in I-DEAS) have this capability. For a more complicated 
geometry Schwarz-Christoffel [7] mapping has been used . The difficulty in 
evaluating integrals involved in the Schwarz-Christoffel transformation 
however makes this technique less attractive. Moreover, mesh generated by 
these techniques may introduce elements with high aspect ratios and 
elements that are highly distorted. 


1.2.2 FREE MESH GENERATION 

This method of generation is best suited for models with 
complicated geometry. SUPERTAB has this capability. The model is 
broken down into sub-areas and sub-volumes. On each of the curves of 
every sub-area and sub-volume the number of elements and their 
concentrations are selected. The software then generates a mesh that is 
consistent with the selected values and satisfies the requirement on the 
aspect ratios and the distortion factors of the elements. 

t NASTRAN MSGMESH is developed by MacNeal-Schwendler Corporation 
□ GIFTS is developed by Sperry Univac Computer System 
• SUPERSAP is developed by Algor Corporation 
1 I-DEAS is developed by Structural Dynamic Research Corporation. 
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Although these pre-processors help in generating acceptable meshes, 
it is still difficult to obtain a mesh that is best suited for the problem at 
hand. The difficulty lies in the definition of the " best " mesh . Is there 
a best mesh for a particular domain ? If so, is there a different one for 
different set of boundary conditions or a different set of loading ? Is 
there a different optimal mesh for different differential equations in the 
same domain ? Answers to these questions are discussed in the following 
sections. 


1.3 OPTIMAL MESH 

Recall from section 1.1 that the functional is a function of the 
unknown or dependent variables. Note that it is also a function of the 
coordinates of the nodes. Therefore it can be expressed as: 

(1.2) 7T = 7r(u 5 ,d k ) 

where, u, is the vector of unknown, d k is the position vector of k th 
node. 


In order to obtain a true minimum on (1.2), in addition to the 
equilibrium equations (1.1), it is necessary to consider the following 
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equations. 


(13) 

ftr 1 SKlj df n 

dd k 2 dd k dd k 


where, r k is the residual vector. 


Solution of (1.3) along with the geometrical constraints will yield the 
optimal locations of the nodes, which when used in (1.1) should result in 
a uj that is closest to u^ct. 


The method seems to be very simple, theoretically. However, the 
non-linear algebraic equations (1.3) are difficult to solve explicitly. Even 
for a simple geometry in one dimension the algebra is very complicated. 
Numerical solutions are also difficult [8]. Some of the solution methods 
for non-linear equations like gradient methods and complex methods have 
been tried but with little success. Among investigators examining this 
problem, Prager[9] has made a note worthy contribution. He examined a 
bar with a linearly varying cross section under tension. He showed that 
the grid producing the desired least potential energy is the one where the 
cross section areas at the nodes form a geometric series. This problem is 
studied in greater detail in the next chapter. 
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1.4 MESH REFINEMENT 


As described in section 1.2 the user needs to select the number of 
nodes and elements in the model. The selection may be the one that 
leads to the best description of the domain geometrically. For example, a 
curved surface could be modelled by a series of interconnected flat 
rectangular facets. The larger the number of facets, the better is the 
model. The selection may also be based upon intuition, past experience 
and engineering judgement. The mesh obtained may be adequate in some 
cases. In other cases, especially when singularities are present, the mesh 
may not be adequate to obtain the results to the accuracy desired. In 
such cases, the meshes need to be refined. 


1.4.1 REFINEMENT PROCESS 


There are three ways of refining a finite element mesh: 

a) The H-method: This method increases the number of elements 
and hence decreases the element size while keeping the polynomial order 
of the shape function constant. 

b) The P-method: This method increases the polynomial order of 
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the interpolation function while keeping the number of elements in the 
model constant. 

c) The R-method: This method redistributes the nodes while keeping 
the element number and the polynomial order of the interpolation 
function constant. 


1.4.1. 1 H - Method 

This method is primarily based on the choice of characteristic length 
of the elements. "Characteristic length " is referred to in a generalized 
sense and is required to define the element topologically. A linear 
element requires one characteristic length, whereas an element of 
rectangular shape requires two characteristic lengths and a triangular 
element requires three characteristic lengths for its definition. In the 
triangular element the three length informations may be any combination 
of lengths and angles. 

Instead of expressing the functional in terms of the position vectors 
of the nodes, as in (1.2), it can be expressed as a function of the 
element characteristic lengths as : 
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( 1 . 4 ) 


7T = ^(Ui.hli) 


where, is the element characteristic length, 1 is the index on the 

characteristic length for element k 

Also, note that there will be geometrical constraints on h^. For 
example, the sum of the element lengths in a particular direction should 
be equal to the overall dimension of the model in that direction. Again 
as described in section 1.3, the function can be minimized with respect to 
the characteristic lengths. 


(1.5) 


dn 

dhik 


1 SKij 

— U; 


2 'dhu, 


u j 


df 

dhik 


Ui = r k = 0 


Solving (1.5) along with the constraints yield the characteristic 
lengths and hence defines the best mesh. Equation (1.5) is equivalent to 
(1.3) cast in the frame work of characteristic lengths. Therefore the 
solution as indicated in section 13 is difficult. A practical procedure 
using this method consists of selecting a coarse initial mesh, solving the 
equilibrium equations and computing the residue r k on each element. The 
set of elements with large values of residues is the region that needs to 
be refined. The identified region can be refined by sub-dividing the 
elements, thus creating new regions, or by deleting all the elements in the 
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region and replacing them by finer elements. However, the new elements 
need to be of the same type as those in the initial mesh. The equations 
of the new model are solved and the residues are computed. If the values 
of the residues are still large, the refinement procedure can be repeated. 
Indeed , it could be used iteratively until the solution meets the 
prescribed accuracy. 

The monotonic convergence of the refinement procedure has been 
studied by Melosh [10] and Key [11]. A convergence theorem has been 
introduced by Carroll and Baker [12], which states : 

Theorem: A necessary consequence of the following refinement sequence 

(1.6) TTji > TTn+l > ' ' • > 7T n+m > ^MOCt 

where, m represents successive refinements of the initial finite element 
mesh n, is the existence of an optimum sub-division such that 

(1.7) 

^n+m ( h ik) < ^n+m (hik) 

where h^ correspond to the optimum mesh. 

The usefulness of this theorem can be explored in the discussion of 
the r - method. The difficulty in using this method is in the estimation of 


10 



the derivatives involved in the computations of the residues. 


1.4.1.2 P - Method 

This method is primarily based on the choice of the order of the 
interpolation function, which in practice, translates to the choice of 
element type. For example, in a two dimensional domain, the basic 
triangular element with three nodes at the three vertices uses a linear 
shape function (p=l). In order to choose quadratic shape functions (p=2), 
the triangular element with six nodes, three at the vertices and three at 
mid-side locations, has to be selected. Similarly, for cubic functions, an 
element with nine nodes is selected. 

Higher order elements generally provide better description of the 
domain geometrically. They are particularly useful in regions where the use 
of lower order elements would result in a mesh with poor aspect ratios in 
those elements. From the point of view of solution accuracy, higher order 
elements are usually more accurate than the lower order elements. But 
this does not mean that increasing the polynomial order indiscriminately 
will always provide point-wise convergence to the exact solution. The 
argument is based on the theory of interpolation. Prenter [13] states that 
this notion on convergence was first dispelled by Meray and later by 
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Runge. He illustrates this with the function f(x) = l/(l+5x2) being 
interpolated by Lagrange polynomial of order 5 and 15 with evenly spaced 
nodes in the interval [-1,1] which display divergence at - 1 and 1. 
Although the example is for a continuous interpolation function rather 
than a piecewise function, as in a finite element model, it shows that 
there is good reason to exercise caution in increasing the polynomial 
order. 

1.4.1J R - Method 

This is a far less explored method. It neither increases the 
polynomial order nor decreases the element character length. The mesh is 
refined simply by re-distributing the nodes in the domain such that the 
potential energy is reduced. 

Recall the theorem introduced by Carroll and Baker, stated in 
section 1.4.1. 1. The theorem indicates that there exists an optimal 
distribution of the nodes in a domain. Any other distribution will yield a 
potential energy higher than the lowest possible for the given number of 
degree of freedom. The theorem also indicates that : given a distribution 
of nodes, a new distribution will be a refinement over the old distribution 
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only if it results in a lower potential energy than the old distribution. 
This fact could be used in an iterative refinement procedures. 


1.4.2 ADAPTIVE MESH REFINEMENT 


The refinement that follows the requirements of the differential 
equation or the boundary conditions closely is called an adaptive 
refinement. This method is used to tailor the mesh, including finer 
elements. It can also provide elements of higher polynomial order where 
necessary as opposed to the method of h or p - refinement all over the 
domain. The practical method mentioned in section 1.4.1. 1 is adaptive. The 
obvious advantage is that it achieves the desired accuracy level while 
keeping the number of degrees of freedom low. 


1.4.3 A - PRIORI AND A - POSTERIORI METHODS 


The classification of methods into a-priori and a-posteriori refers to 
refinement before and after the solution of the equilibrium equations. In 
a finite element program the solution process is one that needs much of 
computer time. If discretization errors can be estimated a-priori, then the 
mesh can be suitably altered to obtain the best accuracy possible by 
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solving the equations only once. Unfortunately there are no practical a* 
priori methods available. The author has not found any in the literature 
survey. This study is an attempt to provide one. There are several a- 
posteriori methods available for refinement. 


1.4.4 USE OF HIERARCHICAL ELEMENTS IN REFINEMENT 


An hierarchical displacement element is one whose stiffness matrix 
contains the stiffness matrices of lower order elements explicitly as 
submatrices[14]. 

Consider a two-node axial element. Its stiffness matrix is given by: 


An hierarchical displacement element of one higher order contains 
an additional node in the middle of the element, as in the conventional 
quadratic element. However the shape function chosen for the midside 
node is different from the one chosen in the conventional element. This 
results in the stiffness matrix: 


A E 
1 


1 -1 

-1 1 

0 0 


0 

0 

16 

3 - 
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Note that the stiffness matrix of the basic element is contained in 
the new matrix as a submatrix. The stiffness matrices of higher order 
elements are built by a similar process if a higher order element is coded 
into the finite element program, it includes stiffness matrices of all lower 
order elements. In the process of refinement if a higher order element is 
chosen, the previously computed stiffness coefficients would still be valid. 
Hence, only a few additional coefficients have to be evaluated. The 
method is easier than the conventional p - method of increasing the 
polynomial order where the computation of the entire higher order 
element stiffness matrix is required. 

Refinement using hierarchical elements is a-posteriori and appears to 
be attractive. However more research work needs to be done in this area. 

1.5 RESEARCH OBJECTIVES 

It is clear that the accuracy of the finite element results is mesh 
dependent. A proper mesh selection procedure is therefore necessary. A 
posteriori methods are adaptive in nature but are expensive in terms of 
computer processing time. On the other hand, a priori methods are not 
adaptive. They use geometrical criterions, element aspect ratios, for 
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example, for improvements. Some of them help estimating the overall error. 
They, however, do not indicate the regions which need refinement. The 
prime objective of this report is, therefore, to develop a criterion which helps 
in identifying the region for refinement or rezoning process even before the 
equilibrium equations are solved. The procedure based on the criterion 
should be able to guide the user in improving the grids. 

Finally, the report itself is based upon the first author’s doctoral 
dissertation [15] at the University of Cincinnati. 



2. ANALYSIS 


2.1 ANALYTICAL APPROACH 


The objective is to develop a practical and efficient procedure of 
grid enhancement and optimization. The thesis is that for many problems 
the minimization of the trace of the stiffness matrix with respect to the 
nodal coordinates, leads to a minimization of the potential energy, and as 
a consequence, provides the optimal grid configuration. To see this, 
consider the governing matrix equation of finite element analysis: 

(2.1) [K]{u> = {f} 

where, [K] is the stiffness matrix, (u) is the array of dependent 
variables, and {f} is the force array. 

Matrix [K] can be viewed as an operator which maps {u} into {f}. In this 
context, since [K] is symmetric, an orthogonal transformation {T}, which 
diagonalizes [K], can be found. That is, 

(2.2) [K] = [T] t [K]|T] 
where [K] is a diagonal matrix. 
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Let [T]{u} and [T]{f} be {u} and {f}. Then the potential energy tt 
may be expressed as: 

(2-3) * - } {u} t (K]{u} - 

In terms of the array components, * becomes: 

( \ 

(2-4) it * t T '‘““I’ ’ fa 

4 = 11 ^ J 

A A 

where the kjj (i=l,2,...,n) are the diagonal entries of [K] 

n ~ „ 

Observe in Equation (2.4) that the last term: £$ u i does not explicitly 

involve the nodal coodinates. Therefore, it does not effect the 
minimization of n with respect to the nodal coordinates. Also, since the 
uf are positive and are independent variables in the minimization of n, 
the minimization of rr with respect to the nodal coordinates occurs when 

A A 

the sum of the ky (the trace of [k]) is a minimum. Since the trace of a 
matrix is invariant under an orthogonal transformation, minimizing the 
trace of [K] is equivalent to minimizing the trace of [K]. 
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2.2 ENERGY APPROACH 


Consider a one - degree of freedom system. The external work done 
(=fu) varies linearly with respect to u. Also the strain energy ( = 1/2 Ku 2 ) 
varies quadratically with respect to u. Potential energy is the difference of 
strain energy and work done. See Figure 2.1. 

From the instant, the structure is loaded the operating point moves 
from the origin to the point where the potential energy reaches its 
minima (equilibrium). 

Now consider the structure with a reduced stiffness (K’< K). The 
new strain energy and the potential energy curves are also shown in the 
figure. Note that the strain energy curve has become slightly flatter. 
Therefore the potential energy curve has reached a new low, which is 
lower than the previous. The displacement has improved from u a to Ub- 
Therefore it is quite clear that in a single degree of freedom case, a less 
stiff structure has a lower potential energy than the stiff one. 

Next, consider an n - degree of freedom system. Using an 
orthogonal transformation matrix, [K] can be diagonalized. This would de- 
couple the degrees of freedom. Therefore each degree of freedom can be 
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compared with the single degree of freedom system as described above. If 

A A 

[K] is the transformed stiffness matrix, then finding minimum k^ 
(i=l,2,...n) would yield the best grid. Since for the mesh configuration the 

-A A -A 

minimum ky have been found, the trace of [K] which is the sum of ky 
will also be a minimum. But the trace is an invariant under orthogonal 
trasformation. Therefore minimization of trace of [K] would lead to 
minimization of potential energy. 

It should be noted that the diagonal dominance of [K] is not 
adversely affected by the minimization of trace. The improved stiffness 
matrix is the result of redistribution of the nodes and of not any arbitrary 
mathematical operation. 
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3. APPLICATIONS 


Applications of concepts developed are illustrated in the examples in 
this chapter. 


3.1 TAPERED BAR 


3.1.1 DESCRIPTION 


Consider the axially loaded tapered bar shown in Figure 3.1. (This is 
the same problem examined by Prager [9] and Masur [16]. The objective is 
to determine a finite element mesh which best predicts the axial 
displacement. Let the bar have length L and let it be divided into n elements 
with n+1 nodes (numbered 0 to n) as shown. Let the areas at the ends of 
the bar be A„ and A,. Let £ be the non dimensional length parameter 
defined as: 

(3.1) x 

e 

L 

Then the area at any particular ( along the bar is 
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(3.2) 


A = Ao(l-ce) 


where c is 
(33) 

Hence, the area at 
(3.4) 

where, & is Xj/L. 

Let the individual elements have a uniform cross section. For 
example, let the i th element have cross section area Aj and length I; as in 
Figure 3.2. (Note that the elements do not necessarily have the same 
length.) Then Aj and 1 5 are: 


c = 


Aq-Ai 

Ao 


the i th node is: 


Ai = A 0 (l-c£i) 


(3-5) 


Ai = 


Aj.i+Aj 

2 


lj = Xj.j = L(£i * ^i-l ) 


The element stiffness matrix for the i t h element is [5] : 
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where E is the elastic modulus. Then the trace r of the global stiffness 
matrix is: 


(3.7) 


T 


n 


2E £ 

*=i 



' Am+V 
.6 ' 


The improved node location occurs when the trace is minimized 
with respect to the nodal coordinates & (i=l,2,...n-l). Hence, by setting the 
derivative of r with respect to & equal to zero we obtain: 

(3.8) 

d t 

Using equation (3.5) and simplifying we obtain 


= 0 = 


^“(Aj-i+Aj) 



6 - 6.1 




(Ai+Aj +1 ) 


Aj.i+Aj 


Ai+A i+1 


6+i* 6 


(6 - 6-i) 2 (6+i- 6) 


(3.9) 


IY l. "l 
ij+1 

2 


f 1. ^ 

*i+l 

2 

ft. 

n+i 


fli+il 

* 

« •> J 

- 1 

+ An 

•i J 

+ cA 0 

L 

k / 


J 

+ 1 


To simplify the analysis it is convenient to introduce the length ratio 



parameter rjj defined as Ij/lj. Then the ratio Ij+j/li may be written as: 


(3.10) 


j w 

iilL _ 11 _ 

1| ” Jl_ Til 

h 


Then L/l, is 


(3.11) 


L . f i . f ft . S;_ 

ii )=i ij j=l r il r il 


n 

where S n is defined as 2 r ji- 

H 


Using this notation Equation (3.9) may be rewritten as 


(3.12) 

Aj+j = A; 


r i+l,l 


til 


+ Aj.x 


r i+l,l 


Til 


cAor^x^ 


r i+l,l 


r il 


+ 1 


Also, & may be written as: 


(3.13) 


£i 


lj + 1 2 + ■ • • + lj 

L 
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ft - 


1 » r 21 + r 81 ♦ • • • + r n 

S» 


Sn 


1 


E 

H 


r ji 



Then, from equation (3.4), A] may be written as 


(3.14) 


Ai 



Specifically, A x and A 2 are 


(3.15) 


Aj = A© 


1 • 



1 - 


(1 * r»i) 

c s n 


To obtain the element area ratios let i=l in equation (3.12). A 2 is 


then 


(3.16) 


A 2 = Ai ^r 2J - 1 j + AqT 2 i + cAor 21 


£ai + I] 

S n 
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Then, by using equations (3.15) to eliminate A x and A 2 , we obtain 


(3.17) 


S»(l ■ f 21 ) = c 


From the first of equations (3.15) we have 


(3.18) 


Aj = Aor 2 i 


or 


Ai_ 

A 0 


= r 21 


Similarly, the second of equations (3.15) leads to 


A 2 = A 0 r 2 i 


(3.19) 


or 



Next, let i=2 in equation (3.12). By using the same procedure we 
obtain 


(3.20) 


r 32 = (1 ‘ r 32 + r 2l) 
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But, in view of equation (3.17) this becomes 


1 - r 32 * (1 - r 21 ) (1 - r 32 + r 21 ) 


( 3 . 21 ) 


or 


IjiOm • r Ji) “ 0 


Thus, 

(3.22) r 32 = r 2 i 


From equation (3.14) we have 


(3.23) 


A 3 


A 0 



A 0 r 21 


Therefore, we obtain 


(3.24) 


A 3 A 2 

a 2 Aj 


Ai 

— — = r 21 = r 32 
A 0 


Proceeding similarly for i=3,4,... we get 
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(3.25) 


A n A n .i 

3 — 

An-i ^u-2 


Thus, we have the relations 


r 31 = r 21 r 32 = r 21 


(3.26) 


3 

*41 = r 43 r 32 r 21 = r 2l 


Til = r 21 


Hence, S; is the geometric series 


(3.27) 


S 5 = 1 + r 21 + r 21 + 


W 1 - r 21 
• + ri>, * — : 


Then, from equations (3.14) and (3.17), we have 


Aj = Aq^I * (1 - r 2 j) J - A 0 r 21 


(3.28) 


and then 


A„ = A 0 r 21 - Aj 
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Then, from equation (33) we see that r 21 is 


(3.29) r„ - (1 • c)'/» 


Finally, by substituting into equation (3.13) we have 


(330) 


Ci = — (1 + r 2 i + r 31 + • • • + r u ) 


(1 - r 3I ) 


f 1 + r 21 + rfj + • ■ ■ + r 2 j ) 


= 1 - Tjl m 1 - (1 - C)*^ 

c c 


This is the result obtained by Prager [9] in his analysis of the same 
problem. 


3.1.2 DISCUSSION 


First, observe that in Equation (330) for a uniform thickness beam 
c is zero and thus $ is undetermined. This means that for a uniform 
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thickness beam the nodal position are arbitrary. That is, all mesh are 
equally optimal for a uniform thickness beam. 

Next, consider again the element stiffness matrix of equation (3.6). 
From equations (3.4) and (3.30) the scalar multiplier is 


(3.31) 


EAj 

T" 


E(Ai.! + Aj) _ A 0 c f l + r 21 ' 
2L(fj - £i-i ) 2L 1 - T 21 


Since this is a constant (independent of i) the element stiffness matrix 
is the same for each element. This means that each element has the same 
strain energy. Masur [15] has suggested that this result is due to the simple 
geometry of the problem. 


Even with this simple geometry however, like the methods discussed 
in section 13, the analysis needed to determine the optimal nodal 
positions has been extremely detailed. With more complex geometries the 
analysis will become intractable. However, it is not necessary to obtain 
recursive relationships by analytical methods as employed in this example. 
The criteria of minimizing the trace of the stiffness matrix is a 
comparatively simple procedure - readily amenable to the development of 
computer algorithms for optimal nodal locations. 
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3.13 NUMERICAL EXAMPLE 


To illustrate the value of optimizing the mesh consider an axially 
loaded bar which tapers to 1/3 the base area as in Figures 3.1 to 3.3. 
Specifically, let P, C, E, and L have the values: 


(3.31) P = 20N, A„ = 0.00 15m J , C = 2/3 
E = 2.07 x 10 u N/m 2 , L = 4m 


The objective is to find the axial displacement. From elementary 
mechanics the axial displacement u at any location x is: 


(3.33) 


u * 


-EL-ln : 

AoEc 


i _£* 

L 


To compare the displacement results of finite element models with 
Equation (3.33), four models of the bar, each having four elements, were 
examined. One of the models had a uniform nodal distribution. Another had 
the "optimal" mesh as developed in Equation (3.33). The remaining two 
models had arbitrarily selected nodal distributions. The nodal 
displacements were evaluated using the four models and compared with 
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the displacements calculated by (333). Table 3.1 shows the results. Table 
3.2 presents an error analysis and also an Lj - norm of the error. As 
expected the optimum mesh produces the least Lj error. 


TABLE 3.1 - Comparison Of Axial Displacements For The Tapered Bar Of 
Figure 3.1 Calculated Using Various Models 


Axial 
location 
x, m 

Exact 

displacement 
eq. (3.30), 
10' 9 m 

Displacements using various models, 10‘ 9 m 

Uniform 

Mesh 

"Optimum" 

Mesh 

(Prager) 

Mesh 3 

Mesh 4 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.5 

33.6276 



33.60639 

33.60639 

1.0 

70.46244 

70.2679 


70.41338 


1.4409860 

106.1461 


105.4820 



2.0 

156.7020 

156.1509 



15.56506 

2.5 

208.3078 



206.8158 

207.1804 

2.5358986 

212.2923 


210.9648 



3.0 

267.8830 

266.5719 




3.3678522 

318.4384 


316.4529 



4.0 

424.5845 

421.1612 

421.940 

417.6195 

417.9814 
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TABLE 3.2 - Error Analysis (Tapered Bar) 


Axial 
location 
x, m 


Error for various meshes, 10' 9 

m 

Uniform 

mesh 

"Optimum" 

mesh 

(Prager) 

Mesh 3 

Mesh 4 

0.0 

0.0 

0.0 

0.0 

0.0 

0.5 



0.02121 

0.02121 

1.0 

0.1945411 


0.04906 


1.4409860 


0.664118 



2.0 

0.5506105 



1.0514 

2.5 



1.492 

1.1274 

23358986 


132769 



3.0 

1311108 




33678522 


1.985439 



4.0 

3.423228 

2.64433 

6.965 

6.6031 

1 2 - norm 

3.7119418 

3.6246009 

7.1232118 

6.780697 


3.1.4 OBSERVATIONS 


[1] The analysis and the numerical results demonstrate the potential 
usefulness of the trace minimization mesh improvement method. 

[2] Minimization of the trace of the stiffness matrix is a relatively 
simple mesh optimization procedure. It is readily adaptable to 
algorithm development. 
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Figure 3.3 - Tapered Bar with Axial Load 





O 
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3.2 HEAT TRANSFER IN AN INFINITE CYLINDER 


3.2.1 CONFIGURATION AND PROBLEM DEFINITION 


Consider an annular cylinder with infinite length having inner and 
outer radii: r 0 and r n . Let the thermal conductivity be k . Let the 
temperatures at the inner and outer radii be: T 0 and T n . Then the 
governing equation for the temperature distribution along a radial lines is: 


(3.34) 


_d_ 

dr 



= 0 


The boundary conditions are: 

(3.35) T = T 0 at r*r 0 

T = T n at r * r n 


The solution of equation 3.34 subject to equation 3.35 is: 


(336) 


T = T n + (T 0 - 


T n ) 
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Next, suppose that the temperature gradient at the inner surface is 


specified as: q 0 . The boundary conditions are then: 

(3.37) ^ = q 0 at r = r 0 

T = T n at r = r n 


In this case the solution of equation 334 is: 


(338) 


T = T n 


r 0 qoln 



l r J 


3.2.2 FINITE ELEMENT FORMULATION AND MESH OPTIMIZATION 


Figure 3.4 shows the finite element model. It consists of a series of 
annular elements. For elements (e) let the inner and outer radii be r e 
and r^. The entries of the stiffness matrix are: 


(3.39) 


kfj = 2 7I7C, / 

r. 


dNf j 


' dNl 

dr 

j 


dr 

k j 
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where are the element conductivity constants and where the 


element shape functions Nf and N® are: 


(3.40) Nf 


( r - • j 

[ r *+i - r «] 


and 



By carrying out the indicated operations the element stiffness matrix 
becomes: 


(3.41) [ kfj ] . S. [1 '/ ] 


where S e is defined as: 


r e + r *-i 

(3.42) S e = 7r«, — 

■ r e-l 


Hence the trace r, of the global stiffness matrix is: 


(3.43) 


r * 2 £ S 

6^1 




where n is the the number of elements. 
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The trace may be minimized with respect to the nodal coordinates 
by setting the partial derivative of r with respect to r e , equal to zero. 
This leads to the relatively simple relation: 


(3-44) 


r«+i _ _r«_ 

r 8 ” r e .x 


By repeated use of this relation the nodal positions are given by: 

r W n 

(3.45) r e = r 0 — 

l r oJ 


3.2.3 NUMERICAL EXAMPLES 

To illustrate the effectiveness of the method considers the annular 
cylinder with the following temperatures specified on the boundaries: 

(3.46) r 0 = 20 mm T 0 = 100° C 

r n = 50 mm T„ = 0° C 
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Let the conductivity be constant throughout the cylinder with value: 1.0. 


Consider two finite elements models, each with four elements: Let 
the first have a uniformly spaced mesh. Let the second have a mesh with 
nodal spacing governed by equation 3.45. The objective is to determine 
the temperature distribution across the thickness. 

The solution of the finite element governing equations lead to the 
results listed in Table 3.3. (The temperatures at the intermediate points, if 
they are not obtained directly, are obtained using linear interpolation 
between the nodal values.) The error is defined as the difference between 
the theoretical results and the finite element results. The mesh governed 
by equation 3.34 (called the "improved" mesh) is found to have zero 
errors at the nodes. Hence, the Lj norm of the errors* is much smaller 
than that of the uniform mesh. 
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TABLE 3.3 - Comparison Between Finite Element Uniform Mesh and Improved 
Mesh Temperature Results with Theoretical Values for 
Temperature Specified Boundary Conditions 


Radius 

(mm) 

Temperature 

C 

Error 

Uniform 

Mesh 

Improved 

Mesh 

Theoretical 

Values 

Uniform 

Mesh 

Improved 

Mesh 

20.0 

100.0 

100.0 

100.0 

0.0 

0.0 

25.148669 

76.21655 

75.0 

75.0 

-1.216559 

0.0 

27.5 

65.354967 

65.920251 

65.24534 

-0.109625 

-0.674909 

31.622777 

50.881148 

50.0 

50.0 

-0.881148 

0.0 

35.0 

39.024743 

39.628662 

38.92596 

-0.098784 

-0.702703 

39.763536 

25.538185 

25.0 

25.0 

-0.538185 

0.0 

42.5 

17.790692 

18.316874 

17.73661 

-0.05408 

-0.580262 

50.0 

0.0 

0.0 

0.0 

0.0 

0.0 

Lj norm of error : 

1.6033656 

1.1340184 


Next, consider the same cylinder but let the temperature gradient be 
specified on the inner boundary. Specifically, let the boundary conditions 
be: 


(3.46) = -5.4567833* C/mm at r 0 = 50 mm 

dr 

T = T n = 0*C at r„ = 50 mm 

(The temperature gradient of -5.4567833 * C/mm on the inner boundary 
leads to the same theoretical temperature distribution as in the first 
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example.) Table 3.4 shows the comparisons between the finite element 
solutions and the theoretical values of the temperature. Once again the 
Lj norm of the error shows that the improved mesh provides better 
results than the uniform mesh. 


TABLE 3.4 * Comparison Between Finite Element Uniform Mesh and Improved 
Mesh Temperature Results with Theoretical Values for 
Temperature / Temperature Gradient Specified Boundary 
Conditions 


Radius 

(mm) 

Temperature 

C 

Error 

Uniform 

Mesh 

Improved 

Mesh 

Theoretical 

Values 

Uniform 

Mesh 

Improved 

Mesh 

20.0 

99.47716 

99.570031 

100.0 

032284 

0.42997 

25.148669 

75.818072 

74.677527 

75.0 

-0.818072 

0322473 

27.5 

65.01327 

65.636818 

65.24534 

0.23207 

-0391478 

31.622777 

50.615125 

49.78502 

50.0 

•0.615125 

0.214981 

35.0 

38.82071 

39.458273 

38.92596 

0.10525 

-0.532313 

39.763536 

25.404668 

24.892511 

25.0 

-0.404668 

0.10749 

42.5 

17.69768 

18.238117 

17.73661 

0.03893 

-0.501507 

50.0 

0.0 

0.0 

0.0 

0.0 

0.0 

Lj norm of error : 

1.245467 

1.0172293 


Since the temperature gradient is specified at the inner boundary, it 
is also of interest to know how the values of the temperature gradients 
obtained using the two finite element meshes compare with each other 


44 











and with the theoretical values. Table 3.5 provides such a comparison, 
(the temperature gradients at the intermediate points, if they are not 
obtained directly, are computed by using forward differences.) The 
improved mesh has a small error at the inner radius as well as a smaller 
norm of errors overall. 


3.2.4 DISCUSSION 


The numerical example show that the "improved" mesh of equation 
3.45 produces results which are closer to the theoretical values than those 
obtained using the uniform mesh. Therefore, the mesh of equation 3.45 is 
an improvement over the uniform mesh for both the temperature fixed 
boundary conditions and for the mixed boundary conditions. 

The values of stiffness matrix traces of the uniform and improved 
meshes are 74.666666 and 70.151974 respectively. This indicates that for 
these examples the trace is not especially sensitive to the nodal locations. 
Therefore, the difference in the Lj norms of error are not great. More 
dramatic difference in the results will occur in problems where the trace is 
more sensitive to changes in the nodal positions. 
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TABLE 3.5 - Comparison Between Finite Element Uniform Mesh and Improved 
Mesh Temperature Gradient Results with Theoretical Values for 
Temperature / Temperature Gradient Specified Boundary 
Conditions 


Radius 

(mm) 

Temperature ’ < 

ri 

Error 

Uniform 

Mesh 

Improved 

Mesh 

Theoretical 

Values 

Uniform 

Mesh 

Improved 

Mesh 

20.0 

-43951853 

-4.8347453 

-5.4567833 

0.861598 

0.622038 

25.148669 

-4.5951853 

-3.8449325 

-4.3396146 

-0.2555571 

-0.4946821 

27.5 

-3.4923413 

-3.8449325 

-3.9685697 

0.4762284 

0.1236372 

31.622777 

•3.4923413 

-3.0577627 

-3.4511702 

-0.0411711 

-0.3934075 

35.0 

-2.816404 

-3.0577627 

-3.1181619 

0.3017579 

0.0603992 

39.763536 

-2.816404 

-2.4317489 

-2.7446192 

-0.0717848 

-0.3128703 

42.5 

-2.3596907 

-2.4317489 

-2.567898 

0.2082073 

0.1361491 

50.0 

-2.3596907 

-2.4317489 

-2.1827133 

-0.1769774 

-0.2490356 

Lj norm of error on gradients : 

1.0986497 

0.9918611 


In summary, the results obtained in this paper confirm that nodal 
positioning by minimizing the stiffness matrix trace leads to either an 
optimum or near-optimum mesh. If the mesh is not optimum it can be a 
starting mesh for other mesh refinement procedures and for procedures 
using element division or element enhancement (h-methods or p-methods). 
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3 3 AIRCRAFT LUG PROBLEM 


33.1 INTRODUCTION 

Accuracy and reliability of finite element computation are among the 
most important considerations in numerical structural analysis. Run time and 
costs are becoming less important. Indeed, the costs incurred in ensuring that 
the results are accurate are negligible as compared with the costs of potential 
consequence of wrong decisions [17]. 

In accurate finite element modelling, a combination of element size 
modification (h-refinement) and element order modification (p-refinement) 
provide the most efficient solution convergence. An exponential rate of 
convergence can be achieved with optimally designed meshes and optimal 
order refinements, [18,19]. From a practical standpoint, however, it is often 
difficult to implement and achieve exponential convergence. Nevertheless, 
Szabo [17] states that "good" results can be obtained through mesh design 
along with element order refinements. 

The most widely used procedure is to: 
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[1] Select large elements where the solution is known to be smooth. 


[2] Select smaller elements where the solution is known to vary rapidly; 
as around points of singularity. 

Szabo suggests that refinement toward the singularity should be in 
geometric progression with the ratio of sizes of about 0.1 to 0.15. These 
are empirical values. When small size ratios are used, bad element aspect 
ratios are often introduced leading to poor mesh design. Extreme element 
aspect ratios are the cause for overestimation of structural stiffness [19]. 

The intention of this study is to find a rationale which will guide 
the analyst in selecting a good mesh. The procedure developed herein has 
shown that an improved mesh can be obtained by minimizing the trace of 
the stiffness matrix. The conventional procedure, as outlined above, is used 
in initial mesh selection. The mesh is then improved by minimizing the 
stiffness matrix trace by moving the grid point locations. The "improved" 
mesh can then be used in the h or p-versi6n of mesh refinement leading 
to a much faster convergence. It is believed that the effort spent in 
minimizing the trace will be rewarded in convergence of the solution after 
fewer h or p iterations. 
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3.3.2 ANALYSIS 


Convergence of finite element solutions are often measured by the 
value of the total strain energy. With improvements in grid, the total strain 
energy usually increases. Mesh improvements based on the total strain 
energy produces convergence in the average sense. An analyst, however, 
generally wants to know the location and the accurate value of the maximum 
stress as well as accurate values of displacement at designated points as 
dictated by design requirements. 

If u is the actual displacement field and if u is the displacement 
predicted by the finite element method, then 1 1 u - u 1 1 is the norm of error 
on displacement. Since the displacement formulation of the finite element 
method renders the structure overstiff, u is larger than u. Therefore u 
converges from below. Hence, from a practical standpoint, a combined 
displacement, stress and strain energy criterion should be used. 

In his study [17] Szabo has used hierarchical basic functions for 
interpolation and has demonstrated p-version convergence in an analysis of 
an aircraft lug as shown in Figure 3.5. The lug is 0.5 inches thick and the rest 
of the dimensions, as shown, are in inches. It is made of isotropic 
material of modulus of elasticity of 30,000 ksi and Poisson’s ratio of 
0.3. The lug is fixed along AB. A circular pin carrying a load of 10.0 kips, 
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inclined at an angle of 45 degrees to the horizontal, fits tightly in the hole of 
the lug and exerts pressure on it. The results of total strain energy changes 
obtained in Szabo’s modellings are listed in Tables 3.6 and 3.7 for easy 
comparison. 

In the present study, the standard QUAD4, QUAD8, TRIA3 and 
TRIA6 elements of MSC - NASTRAN are used. The mesh used by Szabo 
is shown in Figure 3.6. Note that elements 2 and 11 are distorted. In 
attempting to improve the mesh, grid point 8 is moved to a new location as 
shown in Figure 3.7. Next, grid points 2,3,5 and 6 are moved to obtain the 
mesh of Figure 3.8. Table 3.6 shows the effect upon the total strain energy 
along with the reduction in the stiffness matrix trace. (The trace of the 
stiffness matrix is obtained by using a series of DMAP instructions of MSC- 
NASTRAN as indicated in the Appendix.) This improvement is significant 
when compared with results of Szabo (see Tables 3.6 and 3.7) for the same 
number of degrees of freedom. (English units are employed in the tables.) 
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TABLE 3.6 - Comparison of Performance of Meshes 1,2 and 3 



Mesh 1 

Mesh 2 

Mesh 3 

Trace (• 10 a ) 


7.8712 

73158 

Degrees of Freedom 

36 

36 

36 

Strain Energy(* 10* 2 ) 
Displacement at 

1.966856 

2.068929 

2.179273 

node 12 (* 10’ 3 ) 
Von-Mises Stress 

5.52766 

5.92351 

639045 

at node 1 

18.15 

18.5 

1530 

at node 4 

23.42 

24.5 

2031 

at node 18 
MaxJPrin. Stress 

5.663 

5.641 | 

| 

5.618 

at node 1 

-2.056 

-3.446 

-2364 

at node 4 

253 

26.89 

22.24 

at node 18 
Min.Prin. Stress 

4.97 

4.731 

4.696 

at node 1 

-19.09 

-19.98 

-1634 

at node 4 

4.361 

5.803 

4.092 

at node 18 

-1.195 

-1313 

-1328 


The maximum vertical displacement occurs at the tip of the lug (node 
12). In Table 3.6 it is seen that the lower value of trace is associated with a 
higher maximum displacement (node 12). Since finite element models tend 
to present "stiffer than actual" systems, the trace minimization demonstrates 
the mesh improvement. In mesh 2 since only grid point 8 is moved to lower 
the value of trace, the stress values at nodes 1 and 4 are more accurate 
than those predicted by mesh 1. In mesh 3, elements 1 and 3 are 
larger than those of mesh 1. Therefore, their centroids are farther 
away from nodes 1 and 4 which leads to lower stress estimation at 
those nodes. Thus, an improved displacement value is obtained at the 
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expense of stress values. This indicates that the mesh near nodes 1 and 4 
needs to be refined. Mesh 3 can of course be improved by further 
minimization of the trace. 


TABLE 3.7 - Variation of Strain Energy with P-Refinement. 


Polynomial Order 

Degrees of Freedom 

Strain Energy (* 10* 2 ) 

1 

36 

1.73592 

2 

100 

2.76935 

3 

170 

2.83543 

4 

266 

2.85719 

5 

388 

2.86491 

6 

536 

2.86839 

7 

710 

2.86906 

8 

910 

2.86928 


Next, mesh 1 is refined by increasing the polynomial order to 2 to 
obtain mesh la. Mesh 4 is constructed by refining mesh 3 using an h- 
version modification catering to the regions of high stress around nodes 1 
and 4 and around the hole. 

When singularities are not present, an increase in p will result in an 
increase in the rate of convergence. However, if singularities are present, 
p-refinement, with a given mesh, will not necessarily result in an indefinite 
increase in the rate of convergence. However an optimal rate of 
convergence can be obtained by a proper spacing of the mesh around the 
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singularity [26]. Further mesh 4 is improved by reducing the trace of the 
stiffness matrix to obtain mesh 5. Results of the three meshes are listed 
in Table 3.8. 


TABLE 3.8 - Comparison of Performance of Meshes la, 4 and 5. 



Mesh la 

Mesh 4 

Mesh 5 

Trace (* 10 5 ) 

43.166 

32511 

31.402 

Degrees of Freedom 

100 

96 

96 

Strain Energy(* 10' 2 ) 

2.994742 

2.448802 

2.440627 

Displacement at 
node 12 (* 10' 3 ) 

9.02065 

732149 

733578 

Von-Mises Stress 
at node 1 

20.14 

2334 

23.76 

at node 4 

26.79 

30.62 

31.13 

at node 18 

15.22 

22.45 

2035 

Max.Prin. Stress 
at node 1 

-5.086 

- 5.737 

-5.978 

at node 4 

29.46 

33.77 

3438 

at node 18 

16.6 

22.27 

20.54 

Min.Prin. Stress 
at node 1 

-22.19 

-25.68 

-26.18 

at node 4 

6552 

7.805 

8.1 

at node 18 

3326 

-0344 

-0378 


Note that these meshes have almost the same number of degrees of 
freedom and therefore their performances can be compared with each 
other. Since quadratic interpolation is is generally more accurate than 
linear interpolation, mesh la predicts an improved strain energy value 
compared to the other two. However, since the region around nodes 1 
and 4 and also around the hole have small elements in mesh 4 and 5, 
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the stress values predicted are much higher than those of mesh la. Also, the 
displacement values are improved. The trace in mesh 5 is smaller because 
the elements around the hole are less distorted as compared with those of 
mesh 4. This leads to a larger displacement but at the expense of a lower 
stress prediction around the hole. However, since a smaller trace improves 
the mesh in the overall sense, the stress values at 1 and 4 are higher. It 
should also be noted that the strain energy of mesh 5 is smaller than that of 
mesh 4. 

The convergence of finite element results in the energy norm is not 
monotonic [28]. However, a larger displacement indicates a larger work done 
by external loads and consequently a lower potential energy. Again mesh 5 
could be improved by relocating the nodes to lower the trace. 

Next, mesh 5a is constructed by increasing the polynomial order of 
mesh 5 to two. Table 3.9 can be used for a comparative study of the three 
meshes: mesh 1 used by Szabo, mesh la obtained by a p-extension of mesh 
1, and mesh 5a obtained by combined h and p-extension along with the 
improvement procedure based on trace minimization. 

Increases in stress values of mesh 5a from those of mesh 1 vary 
from 40% to 400%. Also, the increase in displacement is 65%. This shows 
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the need for refinements and improvements. A large strain energy value 
of mesh la indicates that a p- version refinement converges faster on an 
average sense. Higher stress values of mesh 5a indicate the superiority of 
the h-version refinement. Moreover, the h-version may introduce distorted 
elements. To obtain the best overall mesh for a given number of degrees 
of freedom, it is therefore useful to improve the mesh by trace 
minimization procedure. 

Conclusions: 


In view of the foregoing results a combined stress, displacement and 
strain energy criterion should be used to monitor the convergence. A 
combined grid improvement and refinement procedure should be used for 
the best results. 

The study confirms that the h and p-extensions lead to improved 
meshes. The study also shows that the stiffness matrix trace is a good 
measure of the quality of a mesh, especially when singularities are present. 
Therefore, the step by step procedure to be followed by an analyst is: 

[1] Select large elements where the solution is known to be smooth. 
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[2] Select smaller elements where the solution is known to vary rapidly. 

[3] Improve the grid by reducing the stiffness matrix trace. 

[4] Perform a numerical solution. 

[5] Refine the grid by using the h-version refinement and by 
introducing smaller elements in high stress areas. 

[6] Improve the grid as in step 3. 

[7] Refine the grid using the p-version refinement. 

[8] Continue to iterate until the convergence criterion is satisfied. 


TABLE 3.9 * Comparison of Performance of Meshes 1, la and 5a. 



Mesh 1 

Mesh la 

Mesh 5a 

Trace (• 10 5 ) 

8.0856 

43.166 

158.43 

Degrees of Freedom 

36 

100 

322 

Strain Energy(* 10' 2 ) 
Displacement at 

1.966856 

2.994742 

3.020364 

node 12 (* 10* 3 ) 
Von-Mises Stress 

5.52766 

9.02065 

9.10035 

at node 1 

18.15 

20.14 

25.25 

at node 4 

23.42 

26.79 

33.71 

at node 18 
Max.Prin. Stress 

5.663 

15.22 

23.48 

at node 1 

-2.056 

-5.086 

-5.897 

at node 4 

25.3 

29.46 

36.91 

at node 18 
Min.Prin. Stress 

4.97 

16.60 

23.90 

at node 1 

-19.09 

-22.19 

- 27.68 

at node 4 

4.361 

6.552 

7.768 

at node 18 

-1.195 

3.326 

- 0.872 
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Figure 3.5- Aircraft Lug 
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Figure 3.6 - Szabo’s Model - Mesh 1 


nj 
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Figure 3.7 - Mesh 2 (Lug Problem) 


nj 



60 



Figure 3.8* Mesh 3 (Lug Problem) 
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Figure 3.9 - Mesh 4 (Lug Problem) 


oj 
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Figure 3.10 - Mesh 5 (Lug Problem) 
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3.4 DISK PROBLEM 


3.4.1 DESCRIPTION 

A disk with a uniform thickness of 0.5 inch and a 20 inch diameter is 
supported at two points B and C on its perimeter. It is loaded at a point A, 
on the perimeter as shown in Figure 3.11. It is modelled using TRIA6 
elements of MSC-NASTRAN. 

As indicated in the previous analysis of the lug, the initial mesh design 
is an important step in the analysis. The circular disk is axisymmetric. If the 
loads are also axisymmetric, then it is advantageous to maintain that 
symmetry by choosing annular ring elements. If the disk is modelled as 
shown in Figure 3.12, then symmetry about four planes is retained. However, 
the boundary conditions and the load warrant only symmetry about one plane 
passing through AD. In this study, the interior nodes are located on the 
circumference of a circle. Let r m be the radius of this circle. The non- 
dimensional parameter £ = r m /r 0 is varied to change the mesh design. 

The graph in Figure 3.13 shows the variation of the trace of the 
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global stiffness matrix and displacement at the center and Figure 3.14 shows 
the strain energy and displacement under the load as a function of £. The 
trace reaches its lowest value at f = £, = 0.53. At values of £ significantly 
different from £„ elements become distorted leading to an increase in the 
trace. As grid points are moved towards the point of application of the load 
and support, the modeling of the region close to the periphery improves. 
This is indicated by the increase in the values of displacement under the load 
and strain energy. However, the improvement is restricted by the increasing 
distortion of the elements which is indicated by the increase in trace values. 
Therefore, both displacement under the load and strain energy reach their 
maximum values at £ = £,. They then decrease in a similar fashion. It is 
evident that £, and £, do not coincide. The point at the center of the disk is 
farthest from the periphery, and therefore, in accordance with the theory of 
Saint - Venant, it should be least affected by the load and the boundary 
conditions imposed. The displacement at the center reaches the maximum 
at £„ the value at which the trace goes to a minimum. This shows that a 
minimum trace procedure yields a good mesh in the regions away from the 
boundaries and loads. However it is not the best mesh for the specific loads 
and boundary conditions applied. In order to achieve this, one has to use h 
and p methods of refinements in the areas close to the boundaries. If the 
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restriction that the interior nodes need to be on the circumference of a 
circle, is relaxed, then the trace minimization procedure would yield a 
mesh that has the least element distorsion. Best results could be obtained 
by iterating on refinement and improvement steps until the error is below 
the tolerance level. 

3.4.2 CONCLUSION 

This study shows that the trace minimization procedure improves the 
mesh in an overall sense. For any specific load and restraint set, other 
mesh refinement techniques need to be used. 
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Figure 3.13 - Graph of Trace and Displacement at the Center of the Disk 
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Figure 3.14 - Graph of Strain Energy and Displacement Under the Load on the 
Disk 
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3.5 LAME PROBLEM 


3.5.1 DESCRIPTION 


One of the classic problems, used as a bench mark by most 
researchers, is that of a cylinder subjected to internal or external pressure. 
Lame provided the theoretical solution for an infinitely long cylinder. Lame’s 
results can be used to evaluate the accuracy of finite element results. Since 
it is possible to have models of only finite length, there is an inherent error 
associated with the model. Moreover, when the cylinder is divided into 
elements, discretization errors are introduced. The focus of this study is on 
the minimization of discretization error by designing good grid patterns. 

Consider a cylinder of inside radius r 0 = 5 cms, outside radius r„ = 10 
cms, and length L = 40 cms as shown in Figure 3.15. Using symmetry of the 
cylinder, only one half of the cylinder is modelled with the nodes on the mid- 
section plane restrained in the axial direction. Two stacks of 
triangular ring elements of MSC-NASTRAN are used as shown in 
Figure 3.16. The nodes common to both stacks are arranged to be 
on a cylindrical surface of radius r m . The non-dimensional parameter 
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£ = r m /r 0 is changed to vary the mesh pattern. 


The graph in Figure 3.17 shows the variation of the trace of the global 
stiffness matrix and the variation of the strain energy with respect to £. 
Figure 3.18 shows that of the average radial displacement at the inside 
surface. The average radial displacement is computed by adding the radial 
displacement values at all the nodes on the inside surface and dividing the 
sum by the number of nodes. The trace reaches its maximum value 
at £ = £, - 0.7072 which is the geometric mean of the inside and the outside 
radii. The strain energy reaches its peak at £ = £,. It is seen that £, is 
slightly smaller than £,. For a uniform mesh £ = £. = 0.745. If strain energy 
is used as a criterion for convergence of the finite element solution, then the 
mesh with £ = £, would provide the best mesh. However, the mesh with 
minimum trace is very close to the best mesh. Moreover it has been obtained 
without solving the equilibrium equations. The study shows that the 
minimum trace mesh is an improvement over the uniform mesh. 

Convergence of strain energy does not guarantee the convergence of 
displacement and stress values [17]. The average radial displacement at the 
inside surface reaches its maximum at £ = £*, where £ d is smaller than £,. 
Once again the mesh with minimum trace is an improvement over the 
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uniform mesh because is closer to £d than £ u . 


3S2 CONCLUSION 

The study confirms that the trace minimization procedure will 
provide a good starting mesh. Fewer refinement iterations will be needed 
to achieve convergence in the finite element solutions. 
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Figure 3.15 * Cylinder (Lame Problem) 
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Figure 3.16- Finite Element Model of the Cylinder (Lame Problem) 
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Figure 3.17- Graph of Trace and Strain Energy (Lame Problem) 
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Figure 3.18 - Graph of Average Radial Displacement at the Inside Surface 
(Lame Problem) 
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4. ALGORITHM DEVELOPMENT 


The example problems have shown that the trace minimization 
procedure yields either an optimal mesh, as in the Prager problem, or a 
near optimal mesh as in the other examples. In either case it yields a 
very good starting mesh. 

In any given problem, it may not be difficult to obtain an expression 
for the trace of the stiffness matrix. It would, however, be very difficult to 
obtain recursive relations by minimizing the trace of the matrix with 
respect to the nodal coordinates in order to obtain the grid configuration. 
Instead, let any arbitrary mesh (usually uniform mesh) be selected. This 
mesh may then be improved by relocating the nodes such that the value 
of the trace is lowered. The algorithm for trace minimization is shown in 
the flow chart in Figure 4.1. 

There are three fundamental issues which are important for the 
success of the algorithm. First, the nodes that should be relocated need 
to be identified. Second, the direction and magnitude of the movement of 
each of the identified nodes need to be determined. Third, a criterion for 
the termination of the improvement iteration loop needs to be established. 
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The node identification step is relatively simple. Let k u and k„ be 
the largest and the smallest diagonal entries on the stiffness matrix. 
Analogous to the methods of bisection, nodes associated with stiffness 
entries larger than 1/2 (ku + k,,) must be relocated. In general, the 

"cut-off" value could be expressed as rj c (Ku + K„), where r? c = 0.5 is one 

specific choice. However, for best results, r? c could be determined by 
numerical experimentation. 

Next, one or several nodes could be moved at a time. The latter 
choice of the two will certainly reduce the CPU time. Let be the set 
of elements connected to node "i". Let be the set of nodes associated 

with elements in <ft. Then the set t/>, can be described as the "neighbor 

set of node i". Note that in FEM, relocation of node "i" will effect a 
change in the stiffness associated with its neighbor set only. One of the 
fundamental requirements for better control in the process is to be able 
to distinguish the effects of each individual change. Therefore two nodes, i 
and j, will qualify for relocation only if there are no common nodes in 
their neighbor sets Vt and 0j • In mathematical terms, the intersection of 
neighbor sets of all qualifying nodes should be an empty set. 

The determination of the direction and magnitude of the movement 
of each identified nodes from its old location to its new location is 
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relatively difficult. Observe that in the Prager problem, the scalar 
coefficient in the determination of the trace (Equation 3.6) is given by 
EA(x)/l(x). Any relocation of the node which increases the element length 
also decreases the cross sectional area. Therefore it reduces the stiffness 
contribution to the trace. Node relocations can be accomplished by 
observing the expressions for element stiffnesses and the role of each of 
the parameters involved such as A(x) and l(x) in the tapered bar problem. 
Each node can be selected and moved to a new location manually by 
using a graphic terminal. However the process is slow, cumbersome and 
inefficient. 

Another approach is, for each identified node, to obtain its neighbor 
set. Then compute the trace of the submatrix corresponding to the 
neighbor set, and store it. The most important step in this approach is 
the determination of the trace gradient. In the one dimensional case, the 
gradient can be computed by difference formulae once the value of the 
trace is known at another point. Therefore, select a new location at some 
distance away. (A discussion on the magnitude of this distance is given in 
the following paragraph). Next, compute the trace of the neighbor set 
corresponding to the new location. The trace gradient computed will 
indicate the direction and magnitude of movement for relocation. In the 
two dimensional case, first compute the trace of the submatrix 
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corresponding to the neighbor set of the node in its original location as 
described before. Second, select a new location for the node in any 
direction and compute the trace. Next, select another location in a 
direction perpendicular to the direction chosen for the selection of the 
first new location. Again compute the trace. Using the three trace values, 
gradients in the two mutually perpendicular direction can be computed, 
which when added vectorially will yield the gradient at the original 
location. Similarly, in the three dimensional case, three new locations on 
three mutually perpendicular lines should be used. Finally, the node 
should be moved in the direction indicated by the gradient. 

The magnitude of movement for gradient computation and for final 
node relocation will be a fixed percentage of the distance between the 
node under consideration and its neighbor in the direction of the 
gradient. However the percentage can be fixed empirically and/or by 
numerical experimentation. 

The objective of this algorithm is to produce a mesh with the least 
trace value. The hypothesis is that the trace minimization procedure will 
distribute the stiffness uniformly among the nodes and elements. Therefore 
the criterion for termination of the improvement iteration loop can be 
based either directly on the decrement in trace value or the uniformity of 
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the stiffness values at all node points. 


To see how uniform stiffness results in a uniform distribution of 
error and therefore yields the best mesh, consider a finite element model 
with n degrees of freedom (d.o.f.). If the mesh is to be refined by 
introducing additional nodes, then it is necessary to know the expected 
improvement in error before a refinement step is undertaken. O.C. 
Zienkiewicz et. al. [20] and Peano et. al. [21] have shown that if the 
n+l th d.oi. is to be introduced hierarchically, then the error in the energy 
norm is: 


(4.1) 



( fn + l 


I^n+l,nUn 

I^n+l,n+l 


) 2 


where, f n+1 is force corresponding to the n+l th d.o.f., K n+lin+1 is the 
stiffness of the n+l th d.o.f., K n+1>n is the off-diagonal stiffiiess relating the 
n+l th d.o.f. to the original n d.o.f. system, and u n is the array of nodal 
displacements of the n d.oi. system. The subscripts n,l of the error e 
refer to the n original d.oi. and the new d.oi. 


Zienkiewicz [22] has used the above error relation to define an 
error indicator in the form: 
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(4.2) 


„2 

*7n,l 


In f N “ +l d J 


K, 


n+l,n+l 


where, f is the finite element residual. 


In an adaptive refinement strategy, these indicators are normally 
calculated for all the d.o.f. corresponding to the next refinement. The 
indicators serve the purpose of identifying the region where refinement is 
necessary. 


Next, the error corresponding to the previous iteration wherein the 
n th d.o.f. was added, is: 


(43) 


Cn-1,1 


2 ( fn ^n,n-l^n-l ) 

= — 


The corresponding error indicator is: 


(4.4) 


Va-1,1 ~ 


/ <r N n dfi 
in 


K, 


n.n 


These derivations are for the hierarchical finite elements. However, 
the error with the conventional finite elements will be of a similar form. 
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The most general method of generating good grids is to have an 
equal distribution of some specified weight function. ( See Eiseman [23] 
for a complete discussion on adaptive grid generation.) Often, the error in 
the finite element solution is used as the weight function [24]. Therefore 
the objective is to distribute the error equally among all elements. 
However, the value of the residual f can be obtained only after the 
equilibrium equations are solved. Nevertheless, one way of obtaining an 
equi-distribution of error a priori is by having uniform element stiffnesses. 
As a consequence, f will be nearly uniform among the elements. The 
trace minimization procedure developed herein produces such a result. 


Consider again the Heat Transfer Example of section 3.2. Note that 
each of the ratios in the optimality condition, Equation (3.44), can be 
equated to a constant 7 . 


(4.5) 


£l _ 

to ’ fi 


r <H-l 

r. 


r n-l 


= 7 


Substituting into Equation (3.42), the element stiffness coefficient is : 
(4.6) S. . 7T k. 

7 - 1 


which is a constant. Therefore the trace minimization procedure 
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produces a uniform element stiffness. 


Finally, observe the graphs of errors in Figures 4.2, 4.3 and 4.4. The 
errors are equally distributed with the improved mesh. There is a skewed 
distribution with the uniform mesh. In order to compare the error 
distribution among the elements, rms errors were calculated using 50 
uniformly spaced points along the length of each element. Also the overall 
rms error for the model was calculated using all the points. Table 4.1 
shows the rms error distribution. Note that in all the cases, the improved 
mesh distributes the error more uniformly than the uniform mesh. The 
rms errors on the elements are almost exactly equal in the case where 
temperatures are specified at the boundaries. Therefore the mesh obtained 
is optimal. Similar results, however, are not obtained in the case where 
both temperature and temperature gradient are specified because of the 
inability of FEM to strongly satisfy the Neumann boundary conditions 
[25]. Nevertheless, it demonstrates the usefulness of the trace minimization 
procedure in a priori grid refinement. 
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TABLE 4.1 - Comparison of RMS Error Distribution Among the Elements 
Between Finite Element Uniform Mesh and Improved Mesh for 
the Two Models with Different Boundary Conditions. 



Temperature B.C. 

Temperature 

/Gradient B.C. 



Temperature 

Gradient 


Uniform 

Improved 

Uniform 

Improved 

Uniform 

Improved 

Element 1 


■ram 

|»| 



Bif ? 7;:* 

Element 2 


1 


1 


■ 

Element 3 

0.4422 

0.5171 


I 



Element 4 

0.2857 

0.5171 





Overall 

■Bl 

0.5210 

0.4494 

0.3586 

0.3755 

0.3781 


Another advantage of using uniform stiffness criterion is that the 
matrix condition number improves. Matrix condition number k may be 
shown that: 


(4.7) 



where is the largest and is the smallest eigen values of matrix K. 
When, definition (4.7) is used, k it is called the spectral condition number. 
If k„ and k u denote the smallest and largest diagonal entry of matrix K, an 
expression for the lower bound for k would be: 

(4.8) k u 

K > 

K. 
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Utku and Melosh [27] have shown that the decimal digits lost in 
computation of displacements in finite elements to be: 

(4.9) D = p + log 10 e q + log 10 /c 

where D is the decimal digits lost in the computation, e q is the residual 
unbalanced error and p is the precision of the machine in decimal digits. 
Note that for a mesh with uniform stiffnesses, k„ « k* and therefore k > 1. 
If K is equal to unity, then digits lost in the computation is the minimum. 
The criterion of uniform stiffness therefore attempts to obtain results with the 
least manipulation error. 
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Figure 4.1 - Flow Chart of Trace Minimization Algorithm 
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Figure 4.1 (cont) 
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Figure 4 .3 - Graph of Error in Temperature in the Cylinder Model with both 
Temperature and Temperature Gradient Specified Boundary 
Conditions 



J0JJ3 



Figure 4.4 - Graph of Error in Temperature Gradients in the Cylinder model 
with both Temperature and Temperature Gradient Specified 
Boundary Conditions 





5. CONCLUSIONS AND RECOMMENDATION 


5.1 CONCLUSIONS 


Making a proper choice of mesh is an important step in the finite 
element analysis for obtaining accurate results. Although engineering 
judgment and prior knowledge of analysis prove helpful in mesh design, 
the procedure of trace minimization makes it possible to obtain good 
meshes without depending upon such judgments and knowledge. Any 
arbitrary mesh may be used. The procedure then relocates the nodes such 
that the value of the stiffness matrix trace is lowered. As a consequence 
it redistributes the total error nearly uniformly among the elements. Thus 
the resulting mesh is either optimal or near optimal. The main advantage 
of the procedure is that a good mesh can be obtained before the 
equilibrium equations are solved. A posteriori methods could be used to 
refine the mesh even further. With the use of the trace minimization 
procedure, fewer a posteriori refinements become necessary to obtain the 
desired accuracy level than when the procedure is not used. 

Most finite element packages have routines to check the correctness 
of the mesh generated. They check node coincidences, element coincidence 


(T - 
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and element distortions. The routines iterate on the elements to check if 
the distortion is lower or higher than a set tolerance. As demonstrated in 
the aircraft lug analysis example, the trace minimization procedure yields a 
mesh with minimum element distortion - a useful byproduct. Moreover, if 
the algorithm presented in Chapter 4 is implemented properly, the CPU 
time required may not be excessively larger than that required for the 
routines used for a distortion check. 

Finally, the procedure improves the stiffness matrix condition number 
and therefore reduces the truncation errors. The benefits derived obviously 
outweigh the cost of implementing the procedure. 


5.2 RECOMMENDATIONS 


The algorithm outlined in Chapter 4 needs to be coded for efficient 
processing. Some of the constants, such as the cut off parameter, need to 
be determined by numerical experimentations. The procedure needs to be 
validated by applying it with problems belonging to classes other than 
those considered in this dissertation work. 

In most structural problems, the stiffness matrix is symmetric and 
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positive definite. The procedure developed is based upon these 
assumptions. It will be useful to develop similar procedures for problems 

with stiffness matrices that are non-symmetric and not positive definite. 
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APPENDIX 


MSC-NASTRAN allows the user to compute intermediate values or 
custom build solution sequences via the Direct Matrix Abstraction 
Program known as DMAP module. These DMAP instructions typically 
preceed the last card "CEND" in the Executive Control Deck. The 
following set of DMAP instructions were used in the trace calculations: 

Nastran Executive Control Deck 

ID PROBLEM 3 , LUG 
TIME 1 
SOL 24 

• 

ALTER 219 

DIAGONAL KGG/KGGD/ $ 

MATGEN , /IDEN/1/NDF/0/0 $ 

DIAGONAL IDEN/IVEC/ $ 

TRNSP IVEC/IVECT/ $ 

SMPYAD IVECT,KGGD, , , ,/TRACE/2////////2 $ 

MATPRN TRACE // $ 

CEND 


Case Control Deck. 
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Where KGG is the global stiffness matrix, KGGD is the vector of 
diagonal stiffness entries, IDEN is an identity matrix of size NDF by 
NDF, IVEC is the vector of diagonal entries of IDEN, IVECT is 
transpose of IVEC and TRACE is a matrix of size 1 by 1 and contains 
the value of the trace. Note that in the instruction for matrix generation 
the variable NDF has to be replaced by the number of degrees of 
freedom in the model. 
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